Thermocapillary effects in driven dewetting and self-assembly of pulsed 

laser-irradiated metallic films 
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In this paper the lubrication-type dynamical model is developed of a molten, pulsed laser- 
irradiated metallic film. The heat transfer problem that incorporates the absorbed heat from a 
single beam or interfering beams is solved analytically. Using this temperature field, we derive the 
3D long- wave evolution PDE for the film height. To get insights into dynamics of dewetting, we 
study the 2D version of the evolution equation by means of a linear stability analysis and by nu- 
merical simulations. The stabilizing and destabilizing effects of various system parameters, such as 
the peak laser beam intensity, the film optical thickness, the refiectivity, the Biot and Marangoni 
numbers, etc. are elucidated. It is observed that the film stability is promoted for such parameters 
variations that increase the heat production in the film. In the numerical simulations the impacts 
of different irradiation modes are investigated. In particular, we obtain that in the interference 
heating mode the spatially periodic irradiation results in a spatially periodic film rupture with the 
same, or nearly equal period. The 2D model qualitatively reproduces the results of the experimental 
observations of a film stability and spatial ordering of a re-solidified nanostructures. 

PACS numbers: 47.54.r, 47.61.-k, 47.55.nb, 81.16.Dn, 81.16.Rf 



I. INTRODUCTION 

Studies of dewetting, rupture and pattern formation in 
thin liquid films are very important for advanced tech- 
nologies and also present an interest for basic physical 
sciences. However, most experimental and theoretical re- 
sults have been obtained for aqueous and polymer films. 

In this paper wc report on our modeling studies of a 
metallic film dewetting by pulsed laser irradiation (PLI). 
In the experiments on laser-irradiated metal films by 
Bischof et alX, both generic scenarios of dewetting and 
rupture were observed, i.e. by the growth of surface per- 
turbations (spinodal regime^ii^ii^) and by nucleation and 
growth of holes^i^. Very recently, several groups reported 
the results of similar experiments on thinner films (thick- 
ness h ~ 3—15 nm) that arc irradiated by a single pulsed 
beamiiSiiSiiiSiii, as in Bischof et ai, as well as by two 
or more interfering pulsed beams (Pulsed Laser Interfer- 
ence Irradiation, PLII)i^ii^ii^. These experiments indi- 
cate that dewetting in such films is primarily spinodal. 

The process of dewetting is analyzed in the cited pa- 
pers using a thermal transport modeling^'^ i^'^i^^ and the 
standard isothermal thin film PDE^ii^. The consen- 
sus is that dewetting is driven by a long-range inter- 
molecular (van der Waals) forces and the thermocapil- 
lary forces, with negligible material evaporation. Most 
interestingly, in the PLII mode the same static interfer- 
ence picture is formed on the film surface at each pulse 
(i.e., the alternating lines of the hot and cold regions in 
two-beams PLII, or the rectangular grid of the hot and 
cold lines in four-beams PLII). Through the thermocap- 
illary fluid flow such lateral spatial nonuniformity of the 
temperature field allows the fabrication of one and two- 
dimensional lattices of metal nanoparticlcs. 

While the general picture has been laid out quite clear 
in the cited papers, there have been no attempts to de- 



velop a consistent PDE-based model of dewetting in a 
metallic film system, which incorporates thermocapillar- 
ity and the spatial and temporal nonuniformitics due to 
laser irradiation. Presentation of such model is the sub- 
ject of this paper. The model allows the studies of dewet- 
ting, rupture and nanopatterning in 3D films for a set 
of laser parameters, such as the pulse shape and repeti- 
tion frequency, the power intensity of the radiation, the 
arrangement and separation distance of the interference 
fringes, and the strength of interference. 

It must be noted here that the influence of radiative 
heating on fluid dynamics was theoretically studied only 
in a handful of paper a^^i^ . Most notably, Oron 
& Pelesi^ studied thermocapillary flows and instabilities 
of evaporative thin aqueous films with constant internal 
heat generation. Oron^S expanded this study to the case 
of irradiation by a single continuous-wave laser beam. 
Grigorie\i2i analyzed mechanisms of passive and active 
feedback control of evaporatively driven instabilities in 
irradiated thin films. The primary conclusion of these 
studies, which are based on long- wave theory, is that irra- 
diation can partially suppress the growth of instabilities. 
Also, Ajaev and Willisi^ studied axisymmetric dewet- 
ting and rupture, in a molten state, of a thin metallic 
film, which has been melted by a single energetic laser 
pulse with a Gaussian spatial shape. They accounted for 
evaporation and long-range intermolecular attraction to 
the substrate and identified the thermocapillary stresses 
as the major driving force of a film evolution. It must 
be noted that neither of these papers considers nonuni- 
form irradiation or the film reflectivity. The latter, as 
has been pointed out in Refsi^ii^ is often the quantity of 
key importance for the dynamics of the ultrathin metallic 
films. 

This paper is organized as follows. In Section |TT] we for- 
mulate the equations of fluid motion for a film irradiated 
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by PLI or PLII, and then obtain the temperature field in 
the long-wave approximation. Using this field, we derive 
the 3D long- wave evolution PDE for the surface height. 
In Section [IIII we perform the linear stability analysis of 
the 2D surface of the film and the numerical simulations 
of the 2D surface dynamics. We show how the surface 
stability, surface shape, rupture time and nanostructurc 
array formation are affected by key dimcnsionless param- 
eters, such as the film optical thickness, peak laser beam 
intensity, number of incident pulses and their duration, 
Marangoni and Biot numbers, and by the lateral spatial 
nonuniformity of the heat production in the PLII mode. 
Also we compare the qualitative and quantitative char- 
acteristics of stability and dewetting in the systems with 
zero and nonzero reflectivity. Section IIVI contains the 
discussion and conclusions. 



II. DERIVATION OF THE TEMPERATURE 
FIELD AND EVOLUTION EQUATION 

In the PLI or PLII experiments, the solid film 
is periodically melted by a laser pulse. Following 
Rcfsi^i^i^i^i^i^iHi^ we disregard periodic cycles of the 
melt resolidification between laser pulses and assume, for 
the purpose of the analysis and nonlinear dynamical sim- 
ulations of the dewetting process, that the metallic film 
is always in the molten (i.e., liquid) state. This is rea- 
sonable, since it has been shown that the mass transport 
in the solid state is negligible compared to the one in the 
liquid state, and so is the film deformation^. In other 
words, our model assumes that between the laser pulses 
the film cools down to the temperature that is just above 
the solidification temperature, and thus it always stays 
liquid. Correspondingly, the dewetting process is treated 
as a continuous one in simulations. 

Thus we assume a thin film of an incompressible New- 
tonian liquid lying on a planar horizontal substrate. The 
mean height of the film, H , is assumed much smaller than 
the lateral dimension i, thus H/L = e ^ 1. The surface 
tension cr is a linear function of the temperature T, 

= 0-„i - 7(f - f,n), f > f„i. (1) 

Here Tm is the melting temperature of the film, dm is 
the surface tension at the melting temperature, and 7 = 
— ~>0. The tildes mark dimensional quantities. The 
governing equations of the melt arc the Navier-Stokes, 
continuity, and the energy equations: 

p{vi + (v • V)v) = V • n + pg, (2) 
V • V = 0, (3) 

pcp (ft + V • Vf ) = kV'T + nj^ + Q, (4) 



where v = {u,v,w) = (mi,{(2,U3) is the velocity field, 
^ij = —p^ij + Tij is the full stress tensor for incompress- 
ible fluid, Tij = p (||i -|- is the viscous stress ten- 
sor, /i is the dynamic viscosity, p is the density, k is the 
thermal conductivity, g is the gravitational acceleration, 
and Q is the internal heat source. The internal heat gen- 
eration is due to the absorption of radiation from the 
monochromatic laser beam. We assume (i) that the free 
surface of the film is optically smooth, non-scattering, 
and non-emissive, and (ii) that the solid substrate sup- 
porting the film is black (which rules out reficctions from 
the film-substrate interface). 

Under the stated assumptions the form of the heat 
source term is given by Bougucr's law22: 

Q = ^^^^ f{i, y, i) exp {S{S ~ h)), (5) 

where / is the laser power intensity, 6 is the spatially 
uniform optical absorption coefficient, h is the position 
of the free surface, R{h) is refiectivity and f{x,y,t) is a 
positive function with mean value one whose functional 
form depends on the interference mode and the temporal 
shape of the laser pulse. 

The boundary conditions at the free surface z — 
h{x, y, t) are: 

(i) The normal and shear stress balances : 

n • 17 • n = -o-V • n + n, (6) 

t-J7-n = t-VCT, ^ ^ \ (7) 

where n is the unit outward normal to the surface, 
t is the unit tangent vector, and 11 is the disjoining 
pressure. The latter is specified in the form 11 = 
{A/6TT)h~^+Bh~^ where A and B are the Hamaker 
constants. The first term is due to the dispersion 
(van der Waals) forces and the second term is due 
to the contributions from the kinetic energy of the 
electrons^^. 

(ii) The kinematic condition: 

w = hf. + uhx -\- vhy, (8) 

that balances the normal component of the liquid 
velocity with the speed of the interface. 

(iii) The temperature boundary condition is given by 
the Newton's law of cooling: 

KTi = -ah {f - fa) , (9) 

where ah is the heat transfer coefficient and is 
the air temperature. 
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The boundary conditions for velocity of the melt flow 
at the solid boundary z = (the substrate) are no-slip, 
u = V ^ 0, and no-penetration, li ~ 0. Two types of the 
temperature boundary condition at the substrate will be 
considered. 

The boundary condition of the first type (TBCl) is the 
Newton's law of cooling at z = 0: 



(10) 



where ag is the heat exchange coefficient and Tg is the 
temperature of the substrate. Note that the usual case of 
a fixed temperature at the solid boundary, T = Ts (which 
corresponds to perfectly conducting substrate) can be 
easily obtained. It suffices to take the limit ctg — > oo 
(or equivalently, the limit /?s — *■ oo, see the definition of 
Ps below) in the dimensionless expressions for the tem- 
perature and in the evolution PDE for the film height. It 
is clear from the form of these equations that this limit 
exists and that taking the limit results simply in some 
terms dropping out of the equations. 

The boundary condition of the second type (TBC2) is 
the continuity of the temperature and the thermal flux 
at z = 0: 



T = 9, kTs = K„ 



(11) 



where 9 is the temperature field in the substrate. The 
substrate is assumed thin, ~ H (where Hg is sub- 
strate thickness). Thus the temperature field 9 can be 
derived in the lubrication approximation. 



A. Nondimensionalization 

We use the following scalings to nondimensionalize the 
problem^!: x = ccL, y = yL, z = zH, h = hH, 
u = uU, i ^vU, w = weU, i = {L/U)t, f = {IH/k)T, 
p = {fiU/eH)p,U = {fMU/eH)U, a = (/.;7/e)a,7 = 
{fiU K / el H)j , where U is the characteristic fiow veloc- 
ity. Typical values of the material parameters are shown 
in Table I. 



B. Temperature Distribution in the Film 

Upon using the scalings, the dimensionless energy 
equation is 

ePe (Tt + uTx + vTy + wT^) = 

T,, + e2 (T,, + Tyy) + {D/2)f{l - R{h)) exp {D{z - h)) 



-C'Br {ul + ul + vl + vl + wl) 



-Br iul+vl) 



e^Br {wl + wl) 



(12) 



where Pe = pCpUH/K is the Peclet number, Br = 
fiU^/HI is the Brinkman number, and D = SH is the 
optical thickness of a film of the uniform thickness H for 
the incident radiation with the mean penetration length 



5 ^ . Using values of the material parameters from Table 
I gives Pe = 0.019 < 1 and Sr = O.IK 1. Thus in the 
leading order the energy equation is 

T,, + {D/2)f{l - R{h)) exp {D{z - h)) = 0, (13) 

and the energy equation in the substrate is: 

e,, + {D/2)fc^^{D{z-h)) = Q. (14) 



The dimensionless boundary conditions for Eqs. ([13 
and (1141) are 



z = h . 



or 



T,=pg{T-Ts) (TBCl) 



T = d, T^^Ve^ (TBC2), 



e = T. 



(15) 



(16) 



(17) 



(18) 



where (3 = ahH/n and Ps = ctsH/n are the Biot num- 
bers, r = Ks/k is the ratio of the thermal conductivity 
of the substrate to one of the film, and hg = Hg/H. Of 
course, the boundary condition (fT8|) is not needed when 
the boundary condition (TBCl) is used. 

Solution of Eq. (fO)) subject to the boundary condi- 
tions (dni) and ([II]), or the solution of Eqs. and ([Ti[) 
subject to the boundary conditions ([T5[) . ([T7|) . and ([T8[) 
gives the temperature field in the film, 

T{x,y,z,t) = (^l-exp(-m) (^l+if 

/V exp iD{z - h)) + /(I - R{h)) +Ts + 
{Ta -Tg+ F{h)f{l - R{h))) (T + z) p, 

(19) 

where K = -1/2D, 

F{h) = -| + exp i-Dh) + -\-K, (20) 

and T = l/Ps (TBCl), or T = hs/T (TBC2). Note that 
Eq. ([19]) results upon the linearization in /3 of the full 
solution of the problem. This is warranted since /3 ^ 1 
(see Table II). 

Substitution z = h \n Eq. (|19p gives the temperature at 
the free surface: 

T*'') = T{x,y,h,t) = Ts-F{h)f{l-R{h)) + 
(T + h) {F{h)f{l - R{h)) + Ta- Ts)p. 

(21) 

In Figs. [1] and [5] we show the typical contour plots 
of the temperature field in the film of the uniform di- 
mensionless height ft, = 1 at a vertical cross-section in 
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the middle of the domain. These figures were obtained 
with (TBCl), R{h) = and (TBC2), R{h) ^ 0, respec- 
tively. (Remarks: The form of the reflectivity function 
R{h) is shown in Eq. below. Unless noted other- 

wise in the text or in a figure caption, the parameter 
values that are used to obtain all figures in the paper 
are taken from Table 11. Also, the dimensionless quan- 
tities are plotted in all Figures, except Figs. [S] and 
[TT] ) In the top row of both Figures, the film surface 
is heated by the spatially and temporally uniform laser 
beam (/ = 1), while in the bottom row it is heated by 
the laser beam with the uniform temporal intensity dis- 
tribution but nonuniform spatial intensity distribution 
f{x,y) = 1 + 0.1(cos(47r(a; - 0.5)) + cos{4Tr{y - 0.5))), 
corresponding to four-beams PLII. If the film is irradi- 
ated uniformly, the temperature remains uniform in any 
horizontal plane in the film (top rows of the Figures), but 
for spatially nonuniform irradiation the temperature in 
the film follows the shape of the intensity distribution as 
shown in the bottom rows of the Figures. 

If the optical thickness D << 1 the radiation passes 
through the film and such films are called optically thin 
or transparent. If Z) >> 1 the radiation penetrates only 
into a very thin boundary layer adjacent to the free sur- 
face of the film and in this case the film is called optically 
thick or opaque. Note that in the dimensionless energy 
equation, Eq. p^ . / characterizes the variation of the 
heat source in a horizontal plane, whereas the function 
a{z) = exp (_D(z — /i)) describes its dependence on the 
depth of the film, z. Since < z < h, a{z) is the increas- 
ing function of the height z, however for the optically 
thin film the difference a{h) — a(0) is a small quantity. In 
other words, the top and the bottom of the film receive 
approximately same energy from the laser beam. As a 
result the temperature difference across the film is small 
(see the right panels of Figs. [l]and[2]). On the other 
hand, the bottom part of the optically thick film receives 
less energy from the laser beam than the top part, which 
makes the temperature difference across the film larger 
(see the left panels of these Figures). 

Before we proceed further, we state the form of the 
effective film reflectivity that we use in the paper. We 
assume the following form after Refsi^iiS: 



D=5 



D=0.1 



R{h) — ro(l — exp(— a^/i)). 



(22) 



where tq and are the material dependent parameters, 
see Tables I and II. 



N 0.5 



N 0.5 




FIG. f : Contour plot of the temperature field in the film at 
y = 0.5 for the case (TBCI) and R{h) = 0. Top row: the 
surface is heated by the uniform (spatially and temporally) 
laser beam. Bottom row: the surface is heated by the tem- 
porally uniform but spatially nonuniform laser beam. Left 
panel: optically thick film. Right panel: optically thin film. 
Notice very small vertical temperature gradient in all cases 
(dT/dz > 0). 



D=5 



N 0.5 




N 0.5' 



D=0.1 
_i_ 



— 1902,7 






-1902.4 — 






1902.2 




-190P 


1902— 



N 0.5 



C. Derivation of the Evolution Equation for the 
Film Height 



FIG. 2: Same as Fig. (TJ but for the case (TBC2) and nonzero 
R{h) given by Eq. (|22p . Note that vertical temperature gra- 
dient (positive) is larger than in Fig. [T] 



In this section we outline the conventional procedure^! 
of the derivation of the evolution equation. 

The dimensionless Navier-Stokes and continuity cqua- 
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tions read: 

eRe{ut + uux + vuy + wuz) = —px + 

[UxX +'U'yy) +u^z, 

eRe{vt + uVx + vvy + wu^) = —py + 

{VXX + Vyy) + Vzz, 

e'Re{wt + uwx + vwy + wwz) = -Pz + 

{WXX + Wyy) + 

t^wzz - G, (23) 
Ux + Vy + Wz = 0, (24) 

where Re = pUH/^ is the Reynolds number (which is 
of order one in magnitude for the case considered here) , 
and G is the gravity number. In the leading order the 
Navier-Stokes and continuity equations read 

Px = Uzz, Py=Vzz, Pz^-G, (25) 
Ux + Vy + Wz = 0. (26) 

Using the long- wave approximation, \hx\, \hy\ << 1 
and the usual assumption of large surface tension^, 
C — e^G~^ (where C is the capillary number) in Eq. 
([6]), the dimensionless normal stress balance condition is 

-p = Cihxx + hyy) + ^ + ^, (27) 

where A and B are the dimensionless Hamaker constants. 
The dimensionless shear stress balance condition reads 



Uz = (Tx + hx(Tz, Vz = C7y + hyUz 



(28) 



where the dimensionless surface tension is given by cr = 

CTm - l{T - T„i). 

Taking the cross-sectional averages of u and v over 
the film height and integrating Eq. (j26p using the no- 
slip condition, we obtain a more convenient form of the 
kinematic boundary condition: 



ht + (hu)^ + {hv)y = 0. 



(29) 



Integrating the first two equations in Eq. (j25[) twice in z 
and applying the no-slip boundary condition yields 



U = {— - hz)px + ZUz\z=h, 

V = {— ~ hz)py + zv^\^=h. (30) 

Cross-sectional averaging of Eqs. ([50)1 results in 
_ 1 /•'^ ^ h ^ 

_ 1 , h , 

^ h J jPy + -^Vz\z=h- (31) 

Since and are expressed in terms of the derivatives 
of the the surface tension, Eq. (pS)) . which in turn is ex- 
pressed in terms of the derivatives of the temperature at 



the free surface, we need to calculate the latter deriva- 
tives. Hence, by differentiating Eq. p9p and evaluating 
the derivatives at z = ft, we find 

tHI) - fii-Rma-Pih+T))F^ih)K^y) 

fR{h)(l-(3{h + T))F{h)hxiy) 
-F{h){l-R{h))fxiy) (32) 
Hh + T)F{h){l^R{h))fx(y)!3, 
Ti''^ ^ [3{f{l-R{h))F{h)+Ta-Ts), (33) 
where 

F,{h) = ]^ + \{TD-l)c^^{-Dh). (34) 

Substituting Eq. (p3|) in Eq. (p8)) through the derivatives 
of cr and using Eq. (P7)) in Eqs. (PT|) . the expressions for 
the average velocities become: 

11 - ^ fix ~h 2 \^xxx H" '^yyx) H" 

A 9B 

(/^ - 3^ - MhTi'^^)hx - MftTi") + 0{Pfx\ 

V — S ^ ^ 3 ['''XXy + h'yyy) "f 

A 2 R 

(■^ - — - MhT^'^^)hy - MhT^'^^ + 0{f3fy), 

(35) 

where Af = 7/2 is the Marangoni number. 

Finally, substituting Eqs. (j35p in the kinematic bound- 
ary condition Eq. and using Eqs. ([55)) results in the 
evolution equation for the film height: 



V 



3 3 \ h 3 



+Mf3{Ta ~ Ts)h'^\'h + MFi{h)f{l - R{h))h^\'h 

+MR'{h)F{h)fh'^\'h 

-Mf3{h + T)R'(h)F{h)fh^\/h 

+ M(3f{\ - R{h)) {F{h) -{h + T)Fi(/i))) h'^Wh] . 

(36) 

Here the prime denotes differentiation. 

In the following sections of the paper, we assume 
the laser irradiation either uniform in the ccy-plane, or 
nonuniform only in the cc-direction. The former situa- 
tion corresponds to the PLI case, i.e. one of a single 
incident laser beam, where the interference is absent. (If 
there is only one beam, we approximate its spatial in- 
tensity distribution on irradiated domain by a constant, 
i.e. / = 1 or / = f{t).) The latter situation corresponds 
to two-beams PLII. In both cases the simpler 2D model 
provides valuable insight into the dynamics of dewetting. 
This model is studied in Section Hill 

It must be noted that in Eq. ([55)) we omitted the term 
-M(l - R{h))V ■ [F{h)h'^\/f] and the similar term pro- 
portional to the small parameter /3. In PLII, when / 
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Physical Parameter 


Typical values 


Dimensionless group 


Definition 


Typical values 


Film thickness (H) 


10 nm 


Scaling parameter (e) 


H/L 


0.01 


Optical absorption coefficient (^) 


10** 


Reynolds number {Rg) 


pUH/fi 


1 


Film density (p ) 


8.92 * 10'' Kgjrr? 


Brinkman number {^Bv^ 




0.11 


T-Tpaf, ppinapitv (r^i 


420 J/KgK 


Peclet number {^Pg^ 


pCpU H f K 


0.019 


Thermal conductivity (film) (/^) 


100 WImK 


Capillary number (C) 


fiU/a 


0.1184 


Thermal diffusivity (^) 


2.675 * 10"^ m^/s 


Gravity parameter (G) 


eoaH^/uU 


3.93 * 10"^^ 


Melting temperature {Tm^ of Co 


1768 K 


Biot number 




10"'' 


Ambient Temperature (Ta^ 


300 7^ 


Biot number (/3s) 


OLsB j K 




Substrate Temperature (T^g) 


1900 Ti' 


Surface tension i^Om) 


f?j^ 1 uXJ 


0.084 


Acceleration of gravity (^) 


9.8 m/s^ 


Marangoni number ( TV/) 


eIH^/2fj.UK 


1.125 * 10"^ 


Dynamic Viscosity (/i) 


4.45 * 10"^ Pa - s 


Hamaker constant [A) 


eA/G-JTiiUH^ 


3.37 * 10"* 


Snrfflpp tptmion (at mpltinp' nt Ii/t^,, 1 


1 88 


Hamaker constant (B) 


eB/fiUH 


1.17 * 10"'^ 




5 * 10"* T/Km'^ 


Melting temperature (Tm) 




1768 


(^Vi arapi'PT'iQfip vplnpifv (TI — u 1 nT~f\ 


PlO 7T7 / Q 


Ambient temperature {To) 




300 


Ppftlf TTii"priQitv" ( 7\ 


10"^ W/m^ 


Substrate temperature (Ts) nTs/IH 


1900 


l~rpfii" trflvmfpr pnpfRpipnf (rvi i 


1 41 * 10* M^/m^/^' 


Optical thickness (D) 


SH 


1 


Heat transfer coefficient (qs) 


» 1 W/m^/C 


Substrate thickness (hs) 


Hs/H 


1 


Hamaker constant (^4) 


1.41 * 10~i** J 


Ratio of thermal 






Hamaker constant {B) 


2.6 * 10"" iV 


conductivities (F) 


f^s j 


1.3 * 10"^ 


Substrate thickness {Hs) 


10 nm 


Pre- factor, 






Thermal conductivity (substrate) (ks' 


1.3 W/mK 


reflectivity function 


ro 


0.44 



Inverse exponent, 
reflectivity function (a^^) 



15.5 nm 



TABLE I: Material parameters. 



TABLE II: Dimensionless parameters. 



III. THE 2D MODEL OF FILM DYNAMICS 



The 2D reduction of the evolution equation ([351) reads: 



is nonuniform in the plane (see Eq. (|43)) ). the omitted 
terms describe the surface shape change by thermocapil- 
larity due to in-planc temperature equilibration by heat 
conduction. In the PLII experiments the in-planc heat 
fluxes are negligible because each pulse last only a few 
nanoseconds, heat losses to the substrate arc large and 
the distance in the plane between the interference fringes 
is much larger than the film thickness. Thus the lat- 
eral temperature profile is approximately static, i.e. it 
is determined by the geometrical arrangement of the in- 
terference fringes. A few sample computations that we 
performed with Eq. ([36]) where the omitted terms are 
present, show that these terms are indeed much smaller 
than the other terms, and their influence on the film dy- 
namics is negligible. This can be also understood by 
noticing that in the experiment (and correspondingly, in 
our modeling) the wavenumber of the surface perturba- 
tion, k, is larger or much larger than the wavenumber, 
q, of the spatially modulated laser light field /. See, for 
example. Fig. [12] which is obtained with fc = 2.2 and 
q ^ 0.157. Thus Vf^qf<s:\7h^ kh, and the terms 
proportional to V/ are smaller than the terms propor- 
tional to V/i. 



ht 



C 



3 



h 3 I 



+M(3{Ta - T,)h'h, + MFi{h)f{l - R{h))h^K 
+MB!(h)F{h)fh^K - Mf3{h + r)R'{h)F{h)fh^K 
+M/?/(l - R{h)) {F{h) -{h + T)Fi{h)) h^h,],. 

(37) 

A. Linear Stability Analysis 

For the purpose of the linear stability analysis, we as- 
sume uniform laser power intensity distribution (/ = 1) 
and a small normal perturbation of the uniform base 
state, h = h^°^ +^{x,t) = 1 + e'^*e*''^, where equals 
to one due to nondimensionalization, and lo represents 
the growth rate of the perturbation having a wave num- 
ber k. Linearizing Eq. p7|) in ^ results in the dispersion 
relation: 



w(fc) 



~ 3^ 



3 

+MR'(1)F{1){-1 + (3{1 
M{l-Ril))i-F,{l)-f3iF{l) 



)e - MI3{Ta - T,)e 
f T))fc2 

-(l + T)Fi(l)))fc2. 

(38) 
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From Eq. ([38)) we find the critical wave number kc such 
that > for < A; < k,.: 



3 f G ^ 



Ts) 



+Mi?'(l)F(l)(-l + /3(l + T)) 
+M(1 - R{1)) (-Fi(l) - 13 (F(l) - (1 + T)F,m)f' 

(39) 

The growth rate w attains its maximum value uimax at 
the wave number km = kc/ V^, which is usually referred 
to as the " most dangerous mode" . 

The terms at the right-hand-side of Eq. ([55]) de- 
scribe: the stabilizing effect of gravity, the stabilizing 
effect of capillary forces, the destabilizing and stabilizing 
effects of the van der Waals component and the elec- 
tronic component of the disjoining pressure, respectively, 
the stabilizing effect of the temperature gradient across 
the film, if Ta > Tg, and destabilizing effect otherwise, 
and the last two terms are due to the volumetric heat 
source. The coefficient of the last two terms in Eqs. 
(IMl) and (IMD, MR'{1)F{1){-1 + (3(1 + T)) + M(l - 
i?(l)) (-Fi(l) - - (1 + T)Fi(l))) is negative for 

the typical parameters values from Table II and for all 
values of D. This holds true for both cases of the bound- 
ary conditions ((TBCl) or (TBC2)) and does not depend 
on the presence of film reflectivity (see Figs. [5]and[6|. 
Thus the term associated with the heat source has a sta- 
bilizing impac t^^'^°i^^ . 



Results for the case (TBCl) and R{h) = 




FIG. 3: Variation of u with k: The dash-dot curve shows uj{k) 
calculated without the term containing the effect of the heat 
source in Eq. (|38|l . the solid curve shows u}{k) calculated with 
all terms included. 
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In this section T = and tq = in Eq. ((22|) . 

The typical graphs of Lu{k) are shown in Fig. [3] The 
solid curve represents w(fc) calculated with all terms at 
the right-hand-side of Eq. (|38p and the dash-dot curve 
shows Lu{k) computed without the last term. (The fifth 
term is automatically zero since R'{1) = due to R{h) = 
0.) As expected both the maximum growth rate and the 
cutoff wave number in the former case are smaller than 
the corresponding quantities in the latter case. 

The Marangoni effect is expressed by the two compo- 
nents, one responsible for the destabilizing thcrmocapil- 
lary effect (the fourth term in Eq. (|38p with Ta < Tg, see 
Table II), and another responsible for the stabilizing ef- 
fect associated with the heat source (the fifth term). Fig. 
U shows the Marangoni effect when the stabilizing com- 
ponent is absent (left panel) and when both components 
are present (right panel). Since the Marangoni number 
M increases with the peak intensity /, the irradiation of 
the film with a high intensity laser has stabilizing influ- 
ence (below the onset of evaporation). 

The impact of the optical thickness D on the stability 
of the perturbation is shown in Fig. [5l It can be seen 
that kc and ujmax both decrease as D increases, while the 
magnitude of the stabilizing effect (i.e., the coefficient) 
increases with D. Thus in agreement with Refsj ^^i^°'^^ , 



FIG. 4: Variation of kc with M: The left panel is the graph 
of kc calculated without the term containing the effect of the 
heat source and the right panel is the graph of kc calculated 
with all terms included. 



as the film becomes more opaque the internal heat gener- 
ation increases and the associated Marangoni effect sta- 
bilizes the film. 

Finally, we discuss impacts of f3s and /?. The cutoff 
wave number increases as Ps increases ( kc ~ 3.1099 
for (3s = 1 and kc ~ 3.2257 for jSg ~ 10"*) approaching 
the constant value ~ 3.23 as f3s — > oo. The maximum 
growth rate is very insignificantly affected by the change 
of (3s , increasing slowly and approaching value 7.65 x 10~^ 
as /3s — > oo. Thus the film is more stable for smaller 
f3s. Similarly, the cutoff wave number kc and the maxi- 
mum growth rate increase insignificantly with increasing 
P (while f3 is kept typically small, as required by Table 
II). Clearly, as /3 or f3s increase the amount of heat in the 
film decreases due to heat losses to the ambient, and the 
film becomes less stable. 



□ 



2 4 6 8 10 

D 



FIG. 5: Variation of kc (left), variation of 
i^max (right), and variation of the coefficient 
M (-Fi(l) - /3 {F{1) - (1 + T)_Fi(l))) (bottom) of the 
last term in Eqs. (|38p and (|39|l with optical thickness D. 
Here T = l/f3s (case (TBCl)). 

2. Results for the case (TBC2) and R{h) / 

In this section T = hs/T and tq = 0.44 in Eq. (|^ . 

The eoefficient of the two stabilizing terms in the 
growth rate equation ([55)1 is plotted in Fig. [HI Unlike 
Fig. [HI the dependence of the coefficient on D is non- 
monotone for both zero and non-zero reflectivity when 
the heat conduction in the thin substrate is taken into ac- 
count. The maximum stabilization occurs in films with 
D 1, and the minimum one in films with D ^ 1 or 

The corresponding maximum growth rate and the cut- 
off wavenumber are shown in Fig. [71 As starts to 
increase from zero, the film becomes more stable. In 
the interval 0.11 < D < 3.45 complete stabilization oc- 
curs, i.e. the growth rate is negative for all wavenum- 
bers. For D larger than 3.45, the growth rate is positive 
for < A: < fee (where kc is shown in the bottom right 
panel), and as D increases the film becomes less stable. 
Such non-monotonous dependence of stability character- 
istics on the optical thickness is, in our opinion, very in- 
teresting and unexpected. When the heat conduction in 
the thin substrate is taken into account, the film heated 
uniformly can be completely stabilized against small per- 
turbations in some interval of the optical thickness pa- 
rameter. Of course, the film is still unstable for pertur- 
bations with amplitudes larger than some critical one. 

Fig. [H] shows the maximum temperature in the film as 
a function of the dimensional film height, for / = 1. (Of 
course, the dimensionless parameters were re-calculated 
for each new value of H.) The temperature is increas- 
ing as the height inreases, confirming the observations 
in Refs J-ii^i^iiiii^ii^iii. The trend in Fig. [5] also signals 
that higher laser power intensity is required to melt thin- 



FIG. 6: Variation of the coefScient, Mi?'(l)F(l)(-l + (3{1 + 
T)) + A/(1-7?(1))(-Fi(l) -/3(F(1) - (1 + T)Fi(l))), of the 
last two terms in Eqs. I|38p and (|39|l with D. Dot curve: 
R{h) = 0; solid curve: R{h) / 0. Here T = hs/T (case 
(TBC2)). The dot curve is included for comparison. 




FIG. 7: Variation of ujmax (top row) and kc (bottom row) 
with D. Dot curve: R{h) = 0; solid curve: R{h) / 0. The 
dot curve is included for comparison. 



ncr filmsiiii^iSiiiiii^ii^iii. Tfic slope of the line is larger in 
the case of R{h) = 0, which shows that neglecting the re- 
flectivity over-predicts the temperature in the film. This 
effect is also discussed in Refi^. 

It can be easily shown that the total heat generated 
in the film, i.e, the integral of the source term Q (Eq. 
([5])) over the film height h, is an increasing function of 
h when 6 is greater than (see also Refs<^ii^). This is 
the reason why the film temperature increases as the film 
thickness increases. 
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FIG. 8: Plot of the maximum dimensional film temperature 
vs. film height. Dot curve: R{h) = 0; solid curve: R(h) ^ 0. 
The dot curve is included for comparison. 

B. Nonlinear Evolution of the Film 

The 2D evolution equation for the film height, Eq. 
([57)) . is solved numerically using the method of lines. 
Integration in time is performed using the ODE solver 
RADAU whereas the discretization in space is car- 
ried out in the conservative form using the positivity- 
preservingS^, second order-accurate finite differencing on 
a spatially uniform staggered grid. 

The perturbation 

h{x, 0) = I + Aocos{kmx), (40) 

where Aq is the small amplitude and km is the wave 
number of the fastest growing instability as identified 
in the linear stability analysis, is imposed as the ini- 
tial condition. Eq. (|37p is solved in the spatial domain 
< a; < 2n7r/fcm (where n is an integer) subject to the 
periodic boundary conditions. 

We performed simulations using different forms of 
f{x,t) that characterizes the spatio-temporal power in- 
tensity distribution of the PLI or PLII at the film surface. 

Firstly, we assume / uniform in both time and space, 
i.e. / = 1. This corresponds to a steady heating of the 
film surface by a single laser beam with a uniform spatio- 
temporal shape. The results of the nonlinear simulations 
in this regime can be compared to the linear stability 
analysis in the previous section. 

Next, with / still spatially uniform we assign a Gaus- 
sian temporal shape to /, so that 

(t-T/2)2 

/ = fit) = e ^ (41) 
for a single pulse laser irradiation, and 

N-l , 

E(t-(2fc + l)T/2)'^ 
e ^ (42) 

for a sequence of A'' pulses of irradiation, where T is the 
pulse duration and a is the standard deviation. (Note 



that d = 2cr\/2 In 2 is the Gaussian full width at half- 
maximum.) Both situations again correspond to heating 
of the film surface by a single laser beam. 

Finally, we consider the case of the two-beam inter- 
ference. The two-beam interference produces a spatially 
modulated light field having the formic 

/ = /(x) = l + acos(q(x--^)), (43) 

where the parameter < a < 1 models the strength of 
the interference and 27r/g = ^ is the distance between two 
neighboring interference fringes. This distance i is given 
by £ = A/2 sin(6'/2), where A is the wavelength of the 
primary laser beam and 6 is the angle between the two 
interfering beams. Note that the subtraction of Tr/fc™ 
from X is to make the beam focused at the center of the 
domain. 

Below we show the numerical results for the case 
(TBCl) and zero reflectivity. The numerical simulations 
for the case (TBC2) and R{h) ^ give similar shape pro- 
files. The inclusion of the reflectivity in the simulation 
makes the rupture time of the film shorter. In this sense 
the film is more unstable when R(h) ^ 0. We discuss this 
issue in more detail at the end of this Section. 

Impacts of the different modes of irradiation 

/ = 1 . Fig. [5] shows the evolution of the film sur- 
face resulting from uniform irradiation, both in time and 
space, by a single laser beam. The minimum point of the 
initial film surface goes further down and touches the 
substrate. In the simulation, the nonlinear dynamics of 
the film is followed until its minimum height gets a very 
small value (0.001) and then the film is said to be rup- 
tured. Favazza et ali^ estimate the rupture time Tr of a 
film (i.e., the time scale of dewetting) as T,. = 27: / uj{km)- 
For a lOnm-thick film their calculation gives Tr ~ 1.25 
ms. Our linear stability analysis gives ~ 1-9 ms. Note 
that Tr depends on the amplitude of the initial surface 
perturbation. Naturally, perturbations with larger initial 
amplitude rupture the film faster. In the nonlinear dy- 
namical simulation shown in Fig. [5] the rupture time for 
^0 = 0.01375 is approximately 0.9 ms, which corresponds 
to the dimensionless time Tr ~ 43150. It is seen in Fig. 
[TO] that at the initial stage of the irradiation the nonlin- 
ear instability growth rate matches the one predicted by 
the linear stability analysis, but towards the end of the 
simulation the nonlinear instability grows much faster. 
This can be explained by observing that when the film 
height gets small the factor A/h in Eq. ([571) becomes 
very large, signaling that the destabilizing van der Waals 
component of the disjoining pressure dominates over sta- 
bilizing forces, which results in a faster instability growth. 
We also noticed that smaller values of the capillary num- 
ber, C, sometimes result in a ring rupture, meaning that 
the film surface touches the substrate some distance from 
the vertical line through the minimum of the perturba- 
tion. Such rupture leaves behind an array of small liquid 
dropsi^. 
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FIG. 9: Profile of the film height (left), and the evolution of 
the minimum point on the film surface (right). The surface 
is heated by a uniform laser beam. 
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FIG. 10: Plots of ln{Ao) + tOmaxt (linear theory, dashed-dot 
curve) and ln{l — hmin) (nonlinear simulation, solid curve) vs 
time. The slope of the solid line equals the growth rate of the 
instability. 



/ as in Eq. The shape profile obtained by irradi- 

ating the surface with a single pulse of width d = 10787 
(200/is dimensional pulse width) is similar to the one re- 
sulting from uniform laser beam heating (Fig. [5]). How- 
ever, Tr w 35138 in this case, which is less than the 
rupture time for the uniform heating. The reason is that 
the energy absorbed from the Gaussian beam is less than 
the energy absorbed from the uniform irradiation, which 
causes the rapid growth of the instability leading to fast 
rupture. 

/ as in Eq. (J^ . 

In this simulation we use several sequences of Gaus- 
sian pulses. In each sequence, pulses have same width 
and also the repetition frequency of pulses is same for 
all sequences. Fig. [TT] shows the rupture time vs pulse 
width. The surface morphology is similar to the one ob- 
tained by uniform laser beam heating. It can be seen that 
the rupture time increases as the pulse width increases. 



Again, this is because the energy generated from a nar- 
rower pulse is smaller than the energy generated from 
a wider pulse, and this increases the growth rate of the 
instability. Heating the film surface with wider pulses 
also require more number of pulses for the film to rup- 
ture. For example, 3064 pulses are enough to get the film 
ruptured when irradiated with a Ins pulse, compared to 
3293 pulses required for a 50ns pulse irradiation. It must 
be noted here that in any irradiation mode the rupture 
time increases approximately linearly with the peak in- 
tensity /. This effect can be seen from Eqs. ([55]) and 
([55)1 . Indeed, larger / does not change the magnitude 
of the fourth (destabilizing) term there (see the defini- 
tions of the parameters M, Ta and Tg in the Table) , but 
it increases the stabilizing contribution of the last term. 
The inclusion in the model of weak evaporation and the 
corresponding recoil pressure at the free surface should 
partially reverse this trend^^. 
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FIG. 11: Variation of the rupture time with the pulse width. 
(The curve is only the guide for the eye.) 

/ as in Eq. |^^ . Fig. [T^] shows the film profile af- 
ter the two-beams irradiation, together with the inter- 
ference field profile, i.e. the function f{x). Values of 
the parameters are a = 0.99, q = 0.157 and k = 2.2. 
The ruptures first occur in the low temperature regions, 
while the strong heat generated in the high temperature 
regions makes the film surface stay near the equilibrium 
position h = I. The rupture time is 600^s, which is com- 
parable to the time estimate given in Ref (> lOOfis for 
a lOnm-thick film). If the irradiation is stopped imme- 
diately after the first ruptures occur, the solidification 
that follows is expected to create the regular array of 
metallic ridges (nanowires) with the axes along the y- 
direction. In contrast to Fig. [9] where the irradiation is 
uniform, here the ridges have different volume, with the 
large (small) volume in the cold (hot) regions. Thus the 
ridges volume distribution has the period £ of the inter- 
ference imprint. This distribution is qualitatively consis- 
tent with the PLII experiments and modelingii. Since 
the bulk film temperature is higher for thicker films (Fig. 
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[8]), then increasing the film height from 10 nm to 15 nm 
ehminates most ridges, except the one in the destructive 
regions (the bottom figure in Fig. [T2|) . (Note that scales 
are very different along the x and h axes in Figs. [5] and 
[T2l thus the true cross-section of each ridge is closer to 
the circular than it appears.) It must be noted that by 
nature of this model the substrate is exposed only at the 
points of film rupture between the ridges in the cold re- 
gions, while in the experiment each ridge terminates at 
the substrate. In other words, the simulation with this 
model can't be continued for t > Tj. and thus we can't 
predict how the mass will be re-distributed if the irradi- 
ation persists. 



ison in the interval D ^ 1. This is because the surface 
is linearly stable for i? = and R ^ Q when D ^ 1, as 
shown in Fig. El and both T^*° and T^=° are formally 
infinite. Thus the ratio is not plotted for 0.085 < D < A.l 
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FIG. 12: Profile of the film height after the two-beams inter- 
ference heating. Top row, left: H = 10 nm, the initial per- 
turbed surface has eight wavelengths in the s-domain; Top 
row, right: H — 10 nm, the initial perturbed surface has 
twenty eight wavelengths; Bottom row: H = 15 nm, the ini- 
tial perturbed surface has twenty eight wavelengths. 

Finally, in Fig. [T3] we plot the ratio of the rupture 
times, T^^°/T^=° vs. D for the case (TBC2). The sim- 
ulations were done, in each case, with the most dangerous 
wavenumbcr km (for R ~ and R ^ the most danger- 
ous wavenumbcrs arc different). It can be seen that this 
ratio is less than one for all values of D and changes 
non- monotonously with D, decreasing first and then in- 
creasing. The minimum value is as small as ~ 10~^. The 
ratio is less than one because the reflectivity reduces the 
heat generation in the film (Eq. thus reducing the 

stabilization. Note that, since the simulation is started 
with the small surface deformation (such that the predic- 
tions of the linear stability theory are valid for at least 
some time), the ratio does not give meaningful compar- 



FIG. 13: Variation of the ratio of the rupture times for the 
simulations with and without reflectivity, with D. (The curve 
in the left panel is only the guide for the eye.) 



IV. DISCUSSION AND CONCLUSIONS 

In this paper we study the dewetting dynamics of a 
pulsed laser-irradiated metallic films. A lubrication-type 
model describing the flow of a molten film and the heat 
conduction in the film is developed. The heat absorbed 
from the laser beam is included as a source term in the 
heat conduction problem and the temperature field dis- 
tribution in the film is obtained by solving this problem 
analytically. In the laser interference mode of irradia- 
tion, we observe that the lateral temperature distribution 
in the film mimics the shape of the lateral power inten- 
sity distribution of the laser. The temperature difference 
across the film and the temperature in the film are higher 
for optically thick films than for the optically thin ones. 

The temperature field is used to derive the 3D long- 
wave evolution PDE for the film height. In order to get 
clear understandings of the film dynamics, we study the 
2D version of the equation by means of the linear stability 
analysis and numerical simulations. 

The linearized problem allows us to investigate the sta- 
bilizing and destabilizing effects of various system pa- 
rameters. Higher peak intensity of the beam and larger 
Marangoni number M either delay the rupture time of 
an initially perturbed film or make the perturbation de- 
cay, while smaller surface Biot number f3 and substrate 
Biot number (3s have the same effect. The increasing 
optical thickness D can have either stabilizing or desta- 
bilizing effect, depending on the magnitudes of the film 
reflectivity and the ratio of the substrate to film thermal 
conductivities. As film becomes thinner, the stabilizing 
effect of the internal heat generation becomes smaller as 
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more heat is generated in thicker films. 

Impacts of the different modes of irradiation are inves- 
tigated numerically in the 2D setting. For the spatially 
uniform (single beam) irradiation the film rupture is spa- 
tially periodic with the wavelength of the fastest growing 
perturbation. The latter wavelength is determined by 
values of all system parameters, including the laser pa- 
rameters. In the two-beam interference heating mode 
the ruptures occur in the (cold) regions of destructive 
interference, while at the (hot) regions of constructive 
interference the initial surface perturbation is still devel- 
oping. Assuming the irradiation is stopped after the first 
ruptures, the solidification is expected to create a ridge 
(nanowire) array, where the spatial distribution of ridges 
volumes follows the spatial periodicity of the interference 
imprint. 

These conclusions do not depend, at large, on the pres- 
ence of the thickness-dependent reflectivity. However, 
quite unexpectedly we found that reflective films with 
D ^ 1 can be completely stabilized against dewetting 
and rupture, although films with D cither small or large 
are less stable than the corresponding non-refiectivc films 
(due to smaller magnitudes of the heat source in the re- 
flective films). The rupture time from the simulations 
is comparable to the estimated and the experimentally 
obtained values^. The rupture time of the films having 
nonzero reflectivity is significantly shorter than the one 
of the non-reflective films. 

Finally, we point out the difference of thcrmocapillary 
mechanisms with and without heat generation in the film 
due to laser beam irradiation of the film surface. In stan- 
dard applications, when the substrate is heated up (i.e., 
the situation where the film surface is not irradiated and 
thus there is no heat production in the film) the thcrmo- 
capillary effect is governed by the M (3{Ta~Ts)h'^V h term 
in the evolution equation, where Ta < Tg. This term is 



responsible for the fluid flow from the high temperature 
region (the one that is closer to the hot substrate, i.e. the 
trough of the surface undulation) to the low temperature 
region, i.e. the crest of the surface undulation, resulting 
in instability and ultimate film break-up in the high T 
region. In that case the temperature gradient across the 
film is negative, dT/dz < 0. On the other hand, when 
the film surface is irradiated and the heat is generated 
in the film, the temperature gradient dT/dz > (see 
Fig. [T|), despite that still Ta < Tg. The heat source term 
in the evolution equation counterbalances the standard 
term, reversing the direction of the fluid flow. In the mul- 
tiwavelength nonlinear simulation shown in Fig. [12] this 
effect manifests as enhanced surface stability in the hot 
regions. The linear stability is also drastically affected 
(sec Figs. 131 in Eland [7]). These results could help under- 
standing the dewetting and rupture process in ultrathin 
metal films irradiated by pulsed laser beams, such as Co, 
Fe, Au, Ni, Cu, Ag and Mo films on glass and Si02/Si 
substratesii'^^^'ii^^iii. 

Future work will focus on development of accurate and 
efficient numerical methods (finite difference and spec- 
tral) for the 3D evolution equation, simulations of the 
corresponding film dynamics, and on quantitative char- 
acterization of the 3D structures size and ordering. 
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